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ABSTRACT 

A  new  wind-input  and  wind-breaking  dissipation  for  phase -averaged  spectral  models  of  wind-generated 
surface  waves  is  presented.  Both  are  based  on  recent  field  observations  in  Lake  George,  New  South  Wales, 
Australia,  at  moderate-to-strong  wind-wave  conditions.  The  respective  parameterizations  are  built  on 
quantitative  measurements  and  incorporate  new  observed  physical  features,  which  until  very  recently  were 
missing  in  source  terms  employed  in  operational  models.  Two  novel  features  of  the  wind-input  source 
function  are  those  that  account  for  the  effects  of  full  airflow  separation  (and  therefore  relative  reduction  of  the 
input  at  strong  wind  forcing)  and  for  nonlinear  behavior  of  this  term.  The  breaking  term  also  incorporates  two 
new  features  evident  from  observational  studies;  the  dissipation  consists  of  two  parts — a  strictly  local  dissi¬ 
pation  term  and  a  cumulative  term — and  there  is  a  threshold  for  wave  breaking,  below  which  no  breaking 
occurs.  Four  variants  of  the  dissipation  term  are  selected  for  evaluation,  with  minimal  calibration  to  each. 
These  four  models  are  evaluated  using  simple  calculations  herein.  Results  are  generally  favorable.  Evaluation 
for  more  complex  situations  will  be  addressed  in  a  forthcoming  paper. 


1.  Introduction 

Over  the  past  four  decades,  the  wave  forecast  is 
routinely  conducted  by  means  of  spectral  numerical 
modeling  of  physical  processes  responsible  for  wave  de¬ 
velopment  and  wave  evolution.  The  evolution  of  the 
wave  spectrum  is  described  by  means  of  the  radiative 
transfer  equation,  which  in  deep  water  can  be  written  as 


dN  X7 
—  +  V- 
dt 


cN 


_  ^ot  -  sin  +  Sds  +  Sn | 


(1) 


where  N  =  N(a ,  6 ,  x,  t)  is  the  wave  action  density  spec¬ 
trum  dependent  on  angular  frequency  a  from  a  frame  of 
reference  relative  to  any  local  currents,  wave  direction  0, 
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distance  vector  x,  and  time  t\  c  is  the  energy  propagation 
velocity  of  the  waves  in  each  dimension;  Stot  represents 
all  energy  fluxes  contributing  to  wind-wave  evolution; 
and  the  action  density  is  related  to  the  energy  density  as 
simply  N=  Eiu.  In  deep  water,  it  is  generally  accepted 
that  wind-wave  growth  is  primarily  a  result  of  three 
physical  processes:  atmospheric  input  from  the  wind  to 
the  waves  Sjn,  wave  dissipation  (resulting  from  breaking 
and  interaction  with  turbulence  and  viscosity)  S^y  and 
nonlinear  energy  transfer  between  the  wave  components 
Snj.  In  finite  depths,  additional  terms  resulting  from  wave- 
bottom  friction  and  triad  interactions  become  significant, 
and  more  terms  can  be  formulated,  which  may  become 
important  in  particular  circumstances.  All  of  these  source 
terms  are  spectral  functions.  The  reader  is  referred  to 
The  WISE  Group  (2007)  for  a  more  complete  overview 
of  this  technique  of  wave  modeling. 

The  input  and  whitecapping  source  terms  used  in  this 
paper  are  based  on  observational  studies  described  by 
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Banner  et  al.  (2000),  Babanin  et  al.  (2001),  Young  and 
Babanin  (2006,  hereafter  YB06),  Donelan  et  al.  (2006, 
hereafter  DBYB),  and  others.  Quantitatively,  the  wind 
input  is  developed  from  direct  measurement  of  surface 
elevation  and  pressure  in  field  conditions,  which  can 
be  expressed  in  terms  of  fractional  energy  increase 
(Donelan  et  al.  2005).  Qualitatively,  there  are  physical 
features  of  behavior  of  both  the  input  and  breaking 
dissipation  observed  in  these  studies,  but  at  present  these 
features  are  still  finding  their  way  into  operational  mod¬ 
eling.  Two  novel  features  of  the  wind-input  source  func¬ 
tion  are  those  that  account  for  the  effects  of  airflow 
separation  (and  therefore  relative  reduction  of  the  in¬ 
put  at  strong  wind  forcing)  and  for  the  nonlinear  be¬ 
havior  of  this  term  (Donelan  et  al.  2005;  DBYB).  The 
whitecapping-dissipation  term  also  incorporates  two 
new  features  evident  from  observational  studies — the 
dissipation  is  a  two-phase  function,  that  is,  it  consists  of 
two  parts,  a  strictly  local  dissipation  term  and  a  cumu¬ 
lative  term  at  smaller  scales  (higher  frequencies),  and 
there  is  a  threshold  for  wave  breaking  below  which 
no  breaking  occurs  (Banner  et  al.  2000;  Babanin  et  al. 
2001;  Babanin  and  Young  2005,  hereafter  BY05;  YB06; 
Babanin,  2009). 

How  arc  these  physics  included  in  the  spectral-model 
source  terms  up  to  date?  In  Wavewatch  III  (Tolman 
2009),  two-phase  behavior  of  the  dissipation  term  is  ac¬ 
commodated,  although  the  assumed  physics  of  Tolman 
and  Chalikov  (1996)  are  different  from  observed  behav¬ 
ior  in  the  experiments  by  BY05  and  YB06,  and  further 
do  not  include  the  threshold  behavior.  Alves  and  Banner 
(2003)  introduced  a  threshold-likc  behavior,  but  the  de¬ 
tails  of  the  implementation  were  not  entirely  consistent 
with  some  known  features  of  the  measurements  (Babanin 
and  van  dcr  Westhuysen  2008).  Both  van  der  Westhuysen 
(2007)  and  Banner  and  Morison  (2010)  implemented 
updated  forms  of  the  Alves  and  Banner  (2003)  dissi¬ 
pation,  incorporating  threshold  limitations  based  on  ob¬ 
servational  values.  However,  these  three  studies  do  not 
apply  two-phase  wave-breaking  dissipation;  the  cumula¬ 
tive  breaking  term  is  omitted,  though  dissipation  simi¬ 
lar  to  that  of  Komen  et  al.  (1984,  hereafter  KHH)  is  used 
to  represent  nonbreaking  dissipation,  that  is,  swell  dissi¬ 
pation.  Thus  far,  the  implementation  of  whitecapping 
dissipation  most  consistent  with  the  qualitative  features 
observed  in  these  studies — Banner  et  al.  (2000),  Babanin 
et  al.  (2001),  YB06,  and  DBYB — was  done  by  Ardhuin 
et  al.  (2010).  Their  formulation  is  different  to  that  of 
B  Y05  and  YB06,  but  they  use  the  same  essential  features: 
a  cumulative  term  is  included  and  the  same  dissipation 
threshold  is  used. 

Surprisingly,  no  essential  attempts  to  implement  the 
DBYB  input  in  a  two-dimensional  (x,  y)  wave  model 


have  been  made  so  far.  For  example,  Ardhuin  et  al.  (2010) 
utilize  the  Janssen  input  (Janssen  1989, 1991),  which  is 
developed  from  the  Miles  theory  (Miles  1957).  Such 
input  does  not  accommodate  the  air-sea  exchange  be¬ 
haviors  observed  by  DBYB  and  Babanin  et  al.  (2007a) 
at  moderate  and  extreme  weather  conditions. 

The  first  implementation  of  the  DBYB  wind  input  and 
the  BY05  and  YB06  whitecapping  dissipation  was  done 
by  Tsagareli  (2009),  and  summarized  in  Tsagareli  et  al. 
(2010)  and  Babanin  et  al.  (2010).  This  was  done  within 
the  WAVETIME  one-dimensional  research  model  with 
the  exact  solution  of  the  nonlinear  integral  in  (1)  (van 
Vledder  2002, 2006).  The  primary  objective  of  the  pres¬ 
ent  paper  is  to  continue  the  work  of  Tsagareli  et  al. 
(2010)  and  Babanin  et  al.  (2010)  by  implementing  the 
new  DBYB  and  BY05  and  YB06  source  functions  in 
a  wave  model  that  can  be  routinely  used  to  produce  two- 
dimensional  forecasts.  Documentation  of  this  stepwise 
development  will  continue  in  future  manuscripts,  most 
notably  with  demonstration  of  practical  applications  of 
the  model. 

The  implementation  of  the  BY05  and  YB06  dissipa¬ 
tion  herein  also  includes  an  improvement  that,  as  will 
be  shown,  allows  much  simpler  calibration,  such  that  it 
does  not  depend  on  wave  age  (discussed  in  section  3). 
The  implementation,  in  addition,  is  generalized  such  that 
the  functional  dependence  of  the  dissipation  on  thresh¬ 
old  exceedence  can,  via  a  parameter  selection,  be  made 
equivalent  to  either  that  used  in  the  other  threshold- 
based  dissipation  modeling  studies  (discussed  above)  or 
that  used  by  BY05  and  YB06.  Thus,  it  is  possible  to 
contrast  in  a  concise  way  the  practical  impact  of  funda¬ 
mental  differences  between  these  previously  proposed 
formulations,  differences  that  until  now  have  not  been 
clarified.  Section  2  describes  the  model  [Simulating  Waves 
Nearshore  (SWAN)],  detailing  the  wind  input  of  DBYB, 
the  dissipation  function  of  BY05  and  YB06,  and  the 
implemented  nonbreaking  dissipation  sink  term,  which  al¬ 
lows  dissipation  to  continue— albeit  much  more  slowly — 
after  the  spectral  density  has  dropped  below  the  breaking 
threshold.  The  new  physics  implementation  is  calibrated 
in  section  3  using  single-point  model  (duration  limited) 
simulations  with  older,  well-known  models  used  as  a 
baseline,  as  opposed  to  standard  fetch-limited  growth 
curve  analysis.  Four  possible  variants  of  the  new  dis¬ 
sipation  source  function  arc  calibrated  in  this  manner. 
The  observation-consistent  dissipation  term  is  then  eval¬ 
uated  in  section  4  by  applying  the  formula  to  parametric 
spectra.  The  behavior  of  the  model  in  free  simulations  is 
evaluated  in  section  5,  first  analyzing  the  time  evolution 
of  the  source  terms,  and  then  a  more  traditional  analysis, 
comparing  wave  growth  against  empirical  expressions. 
Sections  6  and  7  provide  discussion  and  summary. 
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2.  Model  description 

The  third-generation  phase-averaged  wave  model 
SWAN  (Booij  et  al.  1999;  SWAN  Team  2010)  is  se¬ 
lected  as  the  platform  for  this  study,  because,  along  with 
WAVEWATCH  III  (Tolman  1991, 2009),  SWAN  is  one 
of  the  two  such  models  under  active  development  for 
use  by  the  operational  U.S.  Navy.  The  governing  equa¬ 
tion  of  SWAN  is  given  in  (1). 

a.  Wind-input  term 

Tlie  initial  estimate  for  the  wind-input  term  Sin(/,  6) 
is  taken  from  DBYB,  with  three  noteworthy  modifica¬ 
tions.  First,  following  Tsagareli  et  al.  (2010),  a  physical 
constraint  is  applied  to  the  total  stress.  Second,  the  drag 
coefficient  is  changed  to  use  a  recently  proposed  formula. 
Third,  spectral  saturation  is  expressed  in  terms  of  wave- 
number  rather  than  wave  frequency. 

Here  Sin(f,  6)  is  calculated  from  the  following  equa¬ 
tions  (from  DBYB): 


1998),  calculated  from  the  directional  spectrum  nor¬ 
malized  by  the  maximum  value  of  the  spectrum  at  that 
frequency, 


En(m 


Ejm 

£'(/)  ’ 


(8) 


£'(/)  =  max[£(/,0)],0  €  [0,2n],  and  (9) 


A~l 


(/)  =  \2VEn(f,e)de. 

JO 


(10) 


A  constraint  on  the  computed  normal  stress  is  applied, 
which  we  refer  to  as  a  “physical  constraint,”  because  it 
is  a  constraint  on  the  model  physics,  specifically  Sm .  The 
physical  constraint  is  that  the  normal  stress  may  not  ex¬ 
ceed  the  total  stress  less  the  viscous  (tangential  stress) 
and  any  other  stress  that  is  not  estimated,1 


Tnorm  ^tot  ^unknown* 


(ii) 


SinJ(f,e)  =  B(f,  0)£(/,  0),  (2) 

=  ?(/,0)<A  (3) 

P  w 

-K/,0)  =  G  (/,  0),  (4) 

G  =  2.8  -  {1  +  tanh[lO^B„(/)W(/,0)  -  11]},  and 

(5) 


Thus,  the  implemented  constraint  is  Tnorm  <  Ttot  —  tv 
and  we  only  modify  the  wind-input  term  if  Tnorm  >  Ttot  - 
rv.  This  is  done  by  checking  the  initial  estimate  S1IV(/,  9) 
against  the  maximum  possible  normal  stress  level  from 
(11).  For  purposes  of  this  check,  the  nondirectional  form 
of  Sin  is  used,  Sin(f)  =  0)  dO,  following  Tsagareli 

et  al.  (2010).  The  normal  stress  calculation  is 

JlOHzc  (f\ 

~cdf'  (12) 


W(/,0)  =  ^  max 


0,  —  cos(0  -  0  )  -  1 

’  q  V  wv  wn / 


where  /  =  allir  is  the  frequency,  pa  and  pw  are  the  den¬ 
sities  of  air  and  water,  and  6^  and  0wn  are  the  wave  and 
wind  directions,  respectively.  DBYB  calibrate  this  for¬ 
mula  based  on  U  =  U i0,  the  wind  speed  at  10-m  ele¬ 
vation;  thus,  U  ~  U io  is  used  here.  The  parameter  G  is 
the  sheltering  coefficient,  which  accounts  for  the  reduc¬ 
tion  of  atmosphere-to-wave  momentum  transfer  resulting 
from  full  airflow  separation. 

In  a  slight  departure  from  DBYB,  we  use  the  more 
general  form  for  the  spectral  saturation  dimensionless 
variable  in  terms  of  wavenumber,  removing  an  assump¬ 
tion  of  deep  water, 

W)  =  Mf)Bn(f)  =  A(f)(2n)~l *  E(f)k3 *Cg,  (7) 

where  Cg  is  the  group  velocity  and  k  is  the  wavenumber. 
Here  A(f)  is  a  measure  of  narrowness  of  the  directional 
distribution  at  a  frequency  (Babanin  and  Soloviev  1987, 


where  C  is  the  phase  velocity  and  f\  is  the  lowest  modeled 
frequency,  usually  around  0.04  Hz.  The  high-frequency 
tail  of  the  spectrum  is  important  to  stress  calculations, 
so  a  diagnostic  tail  is  added  to  Sjn(/).  Specifically,  Sin(f) 
is  extended  to/ *  10  Hz  using  an  approximation  for  the 
spectral  slope  Sjn(/)  ocf~z  [at  the  high-frequency  limit, 
and  for  E{f)  «  /~5 * *,  the  DBYB  input  has  this  slope  ex¬ 
actly].  The  total  and  viscous  stresses  are  calculated  from 

Ttot  =  U*Pa  =  CDUl0Pa  and  (13) 

Tv  =  CMoPa’  (14> 


1  Values  of  <  60%  are  reported  in  the  literature  as 

being  typical  (Snyder  et  al.  1981;  Hsu  et  al.  1981;  Hsu  et  al.  1982). 

Taken  with  the  fact  that  our  estimated  viscous  stress  ts  generally 

small,  this  suggests  that  these  nonestimated  stresses  may  be  sig¬ 

nificant,  consistent  with  our  presentation  here:  that  Tu,,  —  rv  is  re¬ 

ally  the  upper  limit  on  the  normal  stress  and  may  be  30%-50% 

higher  than  the  actual  normal  stress. 
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where  Cv  is  based  on  the  data  of  Banner  and  Peirson 
(1998),  as  applied  by  Tsagareli  et  al.  (2010), 

Cv  =  max(— 5  X  10"5  £/,„  +  1.1  X  10"3,0).  (15) 

In  this  study,  the  drag  coefficient  CD  is  taken  from 
Hwang  (2011), 

CD  =  1  X  10_4(-0.016f^o  +  0.967 Ul0  +  8.058). 

(16) 


This  formula  is  a  fitting  to  the  datasets  of  Powell  et  al. 
(2003)  and  Jarosz  et  al.  (2007)  and  provides  saturation 
of  the  sea  drag  at  wind  speeds  in  excess  of  30  m  s~~\  and 
even  a  decrease  of  the  drag  at  very  strong  winds.  The 
behavior  is  strictly  empirical,  though  it  has  been  specu¬ 
lated  (Powell  et  al.  2003)  that  one  cause  of  the  reduction 
of  total  stress  could  be  full  airflow  separation,  as  is  ex¬ 
plicitly  included  in  our  wave-supported  stress  imple¬ 
mentation  (5).  Because  there  is  no  upper  limit  on  the 
wind  speed  that  a  user  might  provide  to  SWAN,  it  is 
necessary  to  modify  the  Hwang  formula  to  prevent  the 
drag  coefficient  from  dropping  to  zero  at  extreme  wind 
speeds,  where  reliable  measurements  do  not  exist  to 
provide  guidance.  We  use  the  simplest  possible  modi¬ 
fication,  applied  at  the  t/10,  where  dU*ldU{Q  =  0  in 
Eq.  (16):  U+  =  2.026  m  s”1  for  t/10  ^  50.33  m  s  .  In  the 
cases  presented  herein,  which  are  not  for  hurricane  wind 
speeds,  this  modified  Hwang  expression  yields  drag  co¬ 
efficients  that  do  not  substantially  deviate  from  those 
obtained  with  the  default  formula  of  SWAN,  which  is 
a  modified  form  of  Wu  (1982;  see  section  6  of  WAMDIG 
1988),  though  it  does  necessitate  a  slightly  different  cali¬ 
bration  of  the  dissipation  coefficients  a\  and  a2  (see  sec¬ 
tion  3  here).  This  topic  of  drag  coefficient  will  be  treated 
in  greater  detail  in  a  second  manuscript,  which  does  in¬ 
clude  hurricane  wind  speeds. 

For  the  case  of  rnornl  >  r(ot  -  rv,  it  is  necessary  to 
reduce  the  S*n  to  satisfy  our  constraint.  We  use  the  re- 
duction  Sin_r(f)  =  Lr(f)Sin/f),  where 


M/)  =  min[l, /*(/)]  and 


P(f)  =  exp 


[M)4 


(17) 

(18) 


This  serves  the  same  purpose  as  a  comparable  reduc¬ 
tion  applied  by  Tsagareli  et  al.  (2010).  Here,  the  re¬ 
duction  strength  is  controlled  by  parameter  J?T,  which  is 
determined  via  iteration  such  that  Tnorm  =  rtot  -  rv.  The 
form  of  the  Eq.  (18)  is  designed  to  disproportionately  re¬ 
duce  the  stress  contribution  from  higher  frequencies.  This 
is  based  on  a  practical  consideration — early  experiments 


using  a  model  with  the  DBYB  input  formulation  with¬ 
out  any  reduction  revealed  a  tendency  to  overpredict 
high-frequency  energy.  Also,  as  discussed  in  Tsagareli 
(2009)  and  Tsagareli  et  al.  (2010),  essential  corrections 
to  the  spectral  peak  area  are  to  be  avoided  because  this 
is  where  the  wind  input  was  actually  measured  in  DBYB. 
Thus,  the  implemented  reduction  serves  a  dual  purpose 
of  enforcing  a  physical  constraint  (11),  while  also  re¬ 
ducing  the  overestimation  of  energy  at  high  frequencies. 
An  unintended  consequence  of  the  design  is  that  a  greater 
portion  of  the  reduction  to  Sjn  (/)  is  applied  in  the 
diagnostic  tail  region,  where  S,n  (/)  has  no  influence  on 
model  calculations. 

An  example  of  the  implemented  DBYB  source  func¬ 
tion,  before  and  after  application  of  Lr(/),  is  shown  in 
Fig.  1.  The  default  Sin  (/)  used  by  the  SWAN  model, 
based  on  the  work  of  Snyder  et  al.  (1981)  and  KHH,  is 
also  shown.  The  directional  spectrum  used  as  input  to 
the  two  formulas  is  identical  and  is  taken  from  appli¬ 
cation  of  the  new  model  in  an  unlimited  fetch  scenario 
(section  3)  for  12  h  with  Uxo  =  12  m  s~2,  which  is  also  the 
wind  speed  provided  to  the  two  formulas  here.  The  most 
important  feature  of  the  DBYB  Sin(f)  is  that,  because 
of  the  dependence  (( UIC )  —  l)2  rather  than  (( LHC )  —  1) 
as  of  the  Snyder-Komen  model,  it  tends  to  be  weaker 
close  to  the  peak  (<2 fp  in  the  example  shown)  and 
stronger  at  higher  frequencies  (>2.5/p  in  the  example 
shown).  As  will  be  shown  later,  this  has  a  profound  im¬ 
plication  on  the  required  behavior  of  the  dissipation  term. 
As  expected,  the  impact  of  the  reduction  via  application 
of  Lr(f)  is  primarily  in  the  diagnostic  tail  (for  SWAN,  this 
is  typically  applied  after  1.0  Hz),  but  there  is  still  a  sub¬ 
stantial  reduction  in  the  prognostic  region.  The  contri¬ 
bution  of  each  frequency  bin  to  the  total  stress  is  plotted 
in  the  lower  panel.  The  contribution  is  nearly  constant 
in  the  saturation  region,  with  the  slight  deviation  being 
due  to  the  “-1”  in  (6),  such  that  it  is,  in  fact,  constant  at 
the  high-frequency  limit. 

b.  Whitecapping  dissipation 

The  observation-consistent  whitecapping  term  S^(fy  6) 
used  in  this  study  is  based  on  that  proposed  by  BY05 
and  YB06  having  two  key  features.  The  first  is  that 
waves  do  not  break  unless  the  spectral  density  at  that 
frequency  exceeds  a  threshold  spectral  density  calculated 
from  the  spectral  saturation  spectrum  (see  Banner  et  al. 
2000;  Babanin  et  al.  (2001,  and  section  1).  The  saturation 
spectrum  is  defined  by  Phillips  (1984)  as  B-g~mk9fl 
N(k).  In  terms  of  frequency  spectra,  the  threshold  spec¬ 
tral  density  is  calculated  as 

(19> 
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Fig.  1.  Source  function  calculations  for  UlQ  =  12  ms  x.  (top)  Energy  density  spectrum  used 
for  calculations  (omnidirectional  form,  whereas  calculations  are  made  on  directional  spectrum). 
Frequency  corresponding  to  C  =  Ui0  is  indicated,  as  is  wave-breaking  threshold  as  described  in 
section  2b.  (middle)  Sin(/).  Snyder  et  al.  (1981)/Komen  et  al  (1984)  model  (SWAN  default;  solid 
black  line),  DBYB  prior  to  application  of  L(f)  (blue  line).  DBYB  after  application  of  L(f)y  that 
is,  tw  =  T^un,  (green  line),  (bottom)  Contribution  of  each  frequency  bin  to  the  total  stress  (fre¬ 
quency  grid  is  logarithmic).  The  separation  between  the  prognostic  spectrum  and  diagnostic  tail 
for  a  typical  SWAN  application  (and  all  applications  presented  herein;  black  dashed  line). 


where  Bnx  is  an  empirical  constant  following  the  inves¬ 
tigation  of  wave-breaking  probabilities  by  Babanin  et  al. 
(2007b),  y/B^K  -  0.035.  Once  this  threshold  is  exceeded, 
the  dissipation  depends  critically  upon  the  level  of  ex¬ 
ceedence,  A (/)  =  E(f)  -  Et  (/)•  An  example  threshold 
spectrum  Ej{f)  is  shown  in  the  top  panel  of  Fig.  1. 

The  second  key  feature  of  the  whitecapping  term  is 
that  it  is  two  phase,  because  it  has  been  hypothesized 
to  be  separable  into  two  distinct  mechanisms  (YB06); 
thus,  there  are  two  separate  dissipation  terms, 


V/’0)  =  W0)  +  r2(f,e).  (20) 


The  first  dissipation  mechanism  is  the  inherent  breaking 
component  T\,  which  accounts  for  breaking  resulting 
from  instabilities  of  waves  at  that  frequency, 


w*)  =  *1  A(f)f 


A  (f) 


E(f)l 


E(f,6). 


(21) 


The  second  breaking  component  T2  accounts  for  the 
dissipation  of  waves  induced  by  the  breaking  of  longer 
waves,  for  example,  via  turbulence  created  by  such 
breaking  events.  It  is  not  expected  that  short  breaking 


waves  will  destabilize  longer  waves.  Thus,  the  calculation 
is  based  on  an  integration  of  the  normalized  threshold 
exceedence  from  the  first  prognostic  frequency  up  to  the 
frequency  in  question.  Thus,  T2  is  a  cumulative  term,  not 
local  in  frequency-wavenumber  space, 


T2(f,0)  —  a2 


E(f,  0). 


(22) 


Here  ax  and  a2  are  coefficients  for  calibration,  discussed 
in  detail  below.  The  exceedence  levels  A (/)  are  nor¬ 
malized  by  a  generic  spectral  density  £(/),  which  can  be 
either  E{f)  or  Ej{f).  With  normalization  by  Er(f),  the 
equations  become  similar  to  that  used  by  Ardhuin  et  al. 
(2008),  which  have  since  been  extended  to  directional 
form  by  Ardhuin  et  al.  (2010).  With  normalization  by 
£(/),  the  equation  becomes  similar  to  that  of  BY05 
and  YB06,  Tsagareli  (2009),  and  Babanin  et  al.  (2010).2 


2  It  is  similar,  but  not  identical,  because  in  those  papers,  the 
normalization  by  E(f)  is  outside  the  integral  of  Tr  The  practical 
impact  of  this  difference  will  be  discussed  in  section  3. 
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FlG.  2.  Qualitative  demonstration  of  expected  impact  of  the 
choice  of  normalization  variable  and  choice  of  power  coefficient  L. 
The  vertical  axis  shows  a  dissipation  quantity  that  has  been  nor¬ 
malized  by  the  total  area  under  each  curve  in  the  region  plotted  to 
aid  viewing  and  roughly  approximate  the  effect  of  calibration. 
Here,  only  one  frequency  is  shown,  and  the  ET  at  this  frequency  is 
5  (units  are  unimportant  here),  so  dissipation  is  zero  for  E  <  5. 


The  normalization  allows  us  to  use  the  power  coeffi¬ 
cients  L  and  A/,  which  can  control  how  strongly  the 
dissipation  term  reacts  to  threshold  exceedence.  Figure  2 
illustrates  the  expected  impact  of  the  choice  of  £(/)  = 
ET(f)  versus  £(/ )  and  choice  of  coefficients  L  and  M.  It 
is  clear  that  using  a  power  coefficient  greater  than  unity 
in  combination  with  the  Ardhuin  form  £(/)  —  ET(f) 
implies  that  the  dissipation  will  become  very  large  as  the 
cxcccdcnce  increases,  whereas  the  form  £(/)  =  £(/)  ap¬ 
plied  by  BY05  and  YB06  demonstrates  less  sensitivity. 
This  is  because  [A (f)IE(f)]L  is  always  less  than  unity, 
whereas  for  L  >  1,  [A (f)/ET(f)]L  can  become  large 
quickly  as  the  exceedence  A (/)  is  increased. 

The  directional  narrowness  parameter  A  is  set  to  unity 
in  all  dissipation  calculations.  This  parameter  was  in¬ 
troduced  by  BY05  to  provide  consistency  with  Banner 
et  al.  (2002),  but  more  recent  analysis  by  Babanin  (2009, 
see  discussion  of  Fig.  5.28)  suggests  that  whenever  the 
cumulative  term  is  included,  the  parameter  may  be  un¬ 
necessary.  Further,  it  may  have  a  negative  impact  on  ac¬ 
curacy  at  the  higher  frequencies.  This  is  discussed  in 
more  detail  in  section  6. 

c.  Nonbreaking  dissipation  (swell  dissipation) 

The  new  S formulation  represents  the  rapid  dissipa¬ 
tion  from  breaking  waves.  In  a  decaying  wind-sea  situa¬ 
tion,  once  the  spectral  density  falls  below  the  breaking 
threshold  ET ,  it  can  no  longer  be  dissipated  via  our  5 as 


formula.  It  is  therefore  necessary  to  also  represent  the 
relatively  slow  dissipation  of  nonbreaking  waves,  as 
would  occur  in  this  situation.  For  this  separate  sink  term, 
we  adopt  the  method  of  Ardhuin  et  al.  (2009).  This 
method  uses  one  of  two  calculations,  depending  on 
whether  the  atmospheric  boundary  layer  at  the  sea  sur¬ 
face  is  either  turbulent  or  laminar.  This  is  based  on 
Reynolds  number  Re  =  (4wsip^«orb/vfl),  where  va  is  the 
viscosity  of  the  air,  aOTb  is  half  the  significant  wave  height, 
calculated  for  the  total  energy  (i.e.,  sea  surface  variance) 
as  tforb  =  2 y/EXoV  and  is  the  significant  orbital  ve¬ 

locity  amplitude  at  the  surface,  related  to  the  rms  orbital 
velocity  amplitude  and  rms  orbital  velocity  as  wsig^  = 

V/2«rms^=2wrms- 

For  Re  >  Recr, 

W/.»)  -  (23) 

6  Hw 

For  Re  <  Recr, 

SSwe..(/.0)  -  2^CdSiV*\/2 (24) 

In  simulations  herein,  wc  use  Recr  =  2  X  105,  CdS,v  =  1.2, 
and  fe  of  values  up  to  0.013,  based  on  Ardhuin  et  al. 
(2010).  However,  whereas  those  authors  allow/,,  to  vary 
with  wind  speed  and  Re^  to  vary  with  wave  height,  we 
specify  each  as  constant  for  a  given  simulation. 


3.  Preliminary  model  calibration 

The  wind-input  source  function  Sjn(/)  as  described  in 
section  2a  is  not  calibrated  because  it  is  preferred  that 
the  implemented  source  function  should  not  deviate 
from  the  original  observation-consistent  form  given  by 
DBYB  within  the  frequency  range  /  <  2 fp,  where  mea¬ 
surements  were  made.  The  observation-consistent  S^(f) 
formulation,  on  the  other  hand,  has  four  free  parame¬ 
ters  to  be  set  (alf  a2>  £,  and  A/),  plus  the  selection  of 
£(/)  =  E(f)  versus  £(/) »  ET(f).  The  first  two  pa¬ 
rameters  are  objectively  determined  via  calibration  in 
this  section.  The  latter  three  choices  are  somewhat  more 
subjective,  so  we  present  four  possible  variants  of  the 
model  based  on  these  choices  [L,  A/,  and  £(/)  ],  and 
evaluate  their  merits.  All  calibrations  are  performed 
within  the  SWAN  model. 

Babanin  et  al.  (2010)  adopt  L  —  M  =  1  and  £(/)  = 
£(/).  The  choice  of  L  —  1  is  based  on  observational 
evidence,  namely,  that  T\  should  have  linear  depen¬ 
dence  on  A (/)  (YB06).  Herein,  we  denote  this  form  as 
DL1M1,  where  “D”  indicates  the  concave  down  form 
of  £(/)  =  £(/)  in  Fig.  2.  The  £(/)  =  ET(J)  form  in  Fig.  2 
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is  concave  up  for  all  values  of  L  >  1,  so  this  is  denoted 
with  “U.”  It  is  readily  apparent  from  Fig.  2  that  the 
E(f)  -  ET(f)  model  is  much  more  sensitive  to  L  and 
M  than  the  E(f)  =  E(f)  model.  Thus,  we  include  three 
variants  of  the  former  in  comparisons  herein,  namely, 
UL2M2,  UL1M4,  and  UL4M4,  and  only  one  of  the 
latter,  DL1M1.  For  brevity,  these  four  variants  are  re¬ 
ferred  to  below  as  the  “new  models,”  though  they  are 
all  based  on  a  single  implementation  of  the  new  physics 
in  the  SWAN  code. 

The  traditional  method  of  determining  dissipation 
term  calibration  coefficients  equivalent  to  (a\,  a2)  is  by 
matching  to  fetch-limited  growth  curves.  For  example, 
Tolman  and  Chalikov  (1996)  calibrated  five  free  param¬ 
eters  comparable  to  our  a\%  a2,  L ,  and  M  using  the  Kahma 
and  Calkoen  (1992)  empirical  expressions.  We  diverge 
from  this  tradition  in  two  ways.  First,  we  use  duration- 
limited  wave  growth  rather  than  fetch-limited  wave 
growth.  Only  a  single  grid  point  is  modeled,  zeroing 
propagation  velocities,  in  effect  representing  the  model 
domain  as  an  infinitely  deep  and  wide  ocean  with  no 
spatial  gradients,  V  •  cN  =  0.  This  simple  design  has  an 
obvious  computational  advantage  but  also  makes  it  less 
likely  that  the  results  will  be  contaminated  by  numerical 
error  or  spurious  swell  components  traveling  at  large 
angles  to  the  wind.3  The  second  break  from  tradition 
is  that  we  do  not  use  empirical  growth  curves.  Rather, 
we  compare  the  new  physics  implementation  with  two 
models  that  have  already  been  proven  as  robust  and 
skillful,  at  least  in  the  context  of  wind-sea-dominated 
conditions.  If  results  from  a  new  model  arc  bounded  by 
these  two  thoroughly  tested  models,  then  that  model  can 
be  considered  preliminarily  calibrated,  and  this  further 
constitutes  an  indirect  check  against  historical  observa¬ 
tions  under  nonidealized  conditions.  Further,  results  can 
be  interpreted  in  the  context  of  the  well-known  behavior 
of  these  models.  For  these  two  models,  we  use  the  KHH 
form  as  implemented  in  SWAN,  using  1)  the  dependence 
on  relative  wavenumber  as  favored  by  KHH,  n  —  1,  and 
2)  an  alternate  dependence,  as  used  by  the  author  in 
Rogers  et  al.  (2003)  and  frequently  used  since  then,  n  =  2, 
where  the  KHH  is  S^if)  (klkm)ny  and  km  is  the  mean 
wavenumber.  The  first  KHH  form  is  the  default  setting 
in  SWAN.  The  second  setting  corrects  a  well-known  sys¬ 
tematic  tendency  by  SWAN  (with  the  first  default  set¬ 
ting)  to  underpredict  the  mean  wave  period;  for  further 


3  This  is  a  common  problem  in  idealized  fetch-limited  simula¬ 
tions  using  any  model  without  a  strong  swell  dissipation  mecha¬ 
nism.  Nonlinear  interactions  in  the  model  create  swell  components 
that  are  opposing  the  wind  direction,  that  is,  traveling  toward  the 
zero-fetch  location,  and  the  model  predicts  nonzero  energy  at  this 
location. 


detail,  see  Rogers  et  al.  (2003).  A  time  step  of  30  s  is  used 
in  all  single-point  model  simulations.  The  discrete  in¬ 
teraction  approximation  (DIA)  method  (Hasselmann 
et  al.  1985;  SWAN  2010)  of  computing  Sn\  is  employed; 
implications  of  this  are  discussed  in  section  6.  Testing 
for  realistic,  two-dimensional  applications  is  also  essen¬ 
tial;  this  has  been  done  and  will  be  reported  separately. 

Figure  3  illustrates  the  calibration  process.  [For  math¬ 
ematical  definitions  of  Hm0  and  Tm0 \  and  frequency 
narrowness,  the  reader  is  referred  to  SWAN  (2010).] 
The  case  of  U\0  =  12  m  s-1  is  shown;  other  wind  speeds 
show  comparable  results,  because  the  models  scale  with 
U*  with  only  minor  variations.  Note,  however,  that  be¬ 
cause  the  KHH  models  use  the  modified  Wu  (1982)  Cd 
formulation  (WAMDIG 1988),  they  scale  differently  from 
the  new  models,  most  noticeably  at  high  wind  speeds. 
The  nonbreaking  dissipation  term  Sswc\\(ay  0)  is  not  in¬ 
cluded  in  this  calibration;  this  term  acts  on  the  entire 
spectrum  even  for  wind-sea-dominated  cases,  so  includ¬ 
ing  it  would  introduce  additional  free  parameters  into 
the  calibration  process.  For  this  reason  in  particular,  the 
calibration  herein  should  be  regarded  as  preliminary, 
with  the  final  calibration  to  be  obtained  later  using  com¬ 
parison  to  observations  in  major  ocean  basins  with  both 
wind  sea  and  swell.  Because  5swcn(o-,  0)  is  omitted,  we 
calibrate  the  (a  j,  a2)  coefficients  of  each  model  to  match 
the  higher  wave  height  of  the  two  baseline  models, 
KHH’s  n  =  2.  The  wave  height  comparison  indicates 
one  possible  problem:  energy  growth  during  the  young 
wave  age  stage  of  the  simulation,  with  energy  there 
exceeding  that  of  either  KHH  model.  This  is  attribut¬ 
able  to  the  dependence  of  the  DBYB  wind  input  on 
[(UIC)  -  l]2  and  B'n(o ),  and  is  worst  with  the  DL1M1 
model,  while  it  is  not  noticeable  with  the  UL4M4  model. 
The  7^,01 ,  Tp,  and  frequency  narrowness  comparisons 
mostly  show  good  agreement,  though  it  is  noted  that 
the  DL1M1  model  Tm oi  tracks  closely  the  KHH  n  =  1 
model,  so  we  can  anticipate  that  this  model,  in  real  ap¬ 
plications,  will  tend  to  underpredict  mean  spectral  period 
like  the  KHH  n  —  1  model.  Comparison  of  the  one¬ 
dimensional  spectra  (not  shown)  indicates  a  tendency  of 
these  two  models  to  overpredict  high-frequency  energy. 

These  two  coefficients  (flj,  a2)  can  be  adjusted  in¬ 
dependently,  and  because  a  relatively  strong  T2  will 
generally  increase  the  dissipation  of  higher  frequencies, 
we  find  that  there  is  some  sensitivity  of  mean  spectral 
period  Tm oi  to  the  relative  magnitude  a\\a2.  Thus,  it  may 
be  possible  to  calibrate  the  Ttrto\  growth  curves  of  the 
models  in  this  fashion.  However,  we  do  not  take  this 
approach.  Rather,  we  calibrate  the  a\/a2  ratio  of  each 
model  to  yield  a  frequency-integrated  T\!T2  ratio  that  is 
consistent  with  our  understanding  of  dissipation  in  the 
real  ocean,  namely,  that  for  fully  developed  waves  the 
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FlG.  3.  Calibration  using  singte-gridpoint  simulations.  Calculations  are  for  U\o  a  12  ms  1  and  a  12-h 
duration  is  shown.  Pierson  and  Moskowitz  (1964)  limit  (yellow  dashed  line),  determined  as  by  KHH, 
Booij  et  al.  (1999),  van  der  Westhuysen  et  al.  (2007),  and  others,  with  U*  scaling  = 

5.6  x  10“3;  ££M  =  Ep^fut  =  1 .1  x  103,  with  C*  from  Wu  (1982),  thus  implying  some  deviation  from  the 
value  of  PM,  which  scales  instead  with  U^s.  Frequency  narrowness  is  dimensionless  from  integration 
of  the  one -dimensional  spectrum,  as  defined  by  Battjes  and  van  Vledder  (1984)  and  SWAN  (2010). 


dissipation  should  be  dominated  by  the  induced  term  T2 
(Babanin  et  al.  2010).  This  also  implies  a  certain  level 
of  consistency  between  the  models,  which  aids  in  inter¬ 
preting  later  results.  The  terms  “fully  developed”  and 
“dominated”  are  qualitative  in  nature;  we  semiarbitrarily 
calibrate  (a\,  a2)  such  that  the  contribution  of  T2  after 
12  h  in  the  U\o  =12  m  simulation  is75%-80%  of  the 
total  dissipation.  This  constitutes  a  physical  constraint, 
though  it  is  applied  in  the  calibration  process,  as  op¬ 
posed  to  a  constraint  that  is  enforced  during  every  ex¬ 
ecution  of  the  model,  as  with  the  constraint  on  total 
stress.  Table  1  gives  the  calibration  coefficients  for  each 
of  the  four  models.  Figure  4  shows  the  UL1M4  results 
for  the  first  1.5  h  at  this  wind  speed.  From  the  study  of 
J.-F.  Filipot  (2010,  personal  communication),  which 
utilizes  the  dataset  of  Manasseh  et  al.  (2006),  we  expect 
that  the  transition  at  3 fp  from  majority  T\  to  majority  T2 
dissipation  (i.e.,  the  50/50  crossover  point)  should  hap¬ 
pen  either  at  or  before  the  wave  age  values  correspond 
to  the  Manasseh  et  al.  study  (see  their  Table  1).  Specifi¬ 
cally,  we  expect  the  crossover  to  occur  at  inverse  wave 
ages  U]0iCp  of  3.5  or  larger  (i.e.,  earlier).  The  model  does 


conform  to  this  expectation,  with  the  crossover  occurring 
quite  early,  during  the  first  10  min  of  the  simulation. 

Babanin  et  al.  (2010)  determined  (ai,  a2)  such  that 
the  model  ratio  of  integrated  dissipation  to  integrated 
input  R  =  Dw(IIto(  =  J  S^if)  dfl$  Sin(f)  df  matched  a  pre- 
determined  function  based  on  observations  (Donelan 
1998;  see  section  5a  below).  An  unusual  feature  of  their 
dissipation  is  that  the  (au  a2)  coefficients  were  a  func¬ 
tion  of  wave  age,  varying  by  a  few  orders  of  magnitude 
over  the  typical  wave  age  range.  In  our  own  simulations 
with  the  original  BY05  and  YB06  form  of  the  dissipa¬ 
tion,  we  also  found  that  (au  a2)  must  vary  with  wave  age 
in  order  to  achieve  satisfactory  results,  and  initially 


Table  1.  Coefficients  used  for  each  of  the  four  models.  E  is  the 
normalization  variable,  E(f)  -  E(f)  or  E(f)  =  ET(f). 


E(f) 

L 

M 

0i 

«2 

E(f) 

1 

1 

2.0  x  10'4 

1.6  x  10-3 

EAf) 

2 

2 

8.8  x  10'6 

1.1  x  10~4 

EAf ) 

1 

4 

5.7  x  10-5 

3.2  x  10-6 

EAf) 

4 

4 

5.7  x  10-7 

8.0  x  10~6 
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U  1Q  =  12  nVs  ;  T  and  T2  integrated  over  0.04  to  1 .0  Hz 


00:15  00:30  00:45  01:00  01:15  01:30 

time  (HH:MM) 


Fig.  4.  (top)  Contribution  of  T\  vs  T2  to  the  total  dissipation  for  the  UL1M4  Ui0  =  12  m  s'1 
simulation.  The  integrated  T\  and  T2  values  (solid  lines)  and  the  values  at  3 fp  (dashed  lines), 
(bottom)  Inverse  wave  age  for  the  same  simulation.  Wave  age  and  time  period  corresponding 
to  the  Manasseh  et  al.  (2006)  study  are  indicated  (dashed-line  box). 


developed  complicated  lookup  tables  to  be  used  by  the 
model  to  determine  the  coefficients  during  run  time. 
However,  by  modifying  the  T\(f)  calculation  such  that 
the  normalization  occurs  within  the  integral  as  shown 
in  (22),  we  find  that  a  single  pair  of  coefficients  can  be 
used,  and  the  lookup  table  procedure  becomes  obsolete. 

Summarizing  the  calibration  procedure,  we  find  for 
each  of  four  models  two  coefficients  (aj,  02)  that  yield 
the  best  match  of  total  energy  in  idealized  duration- 
limited  simulations  to  two  preexisting  models,  one  being 
the  default  physics  of  SWAN  (version  40.72).  The  rela¬ 
tive  size  of  ax  versus  ^2  is  set  to  achieve  a  desired  ratio 
of  Tx  versus  Ti  at  large  wave  age  values  at  moderate 
wind  speeds,  thereby  ensuring  some  degree  of  consis¬ 
tency  between  models  with  regard  to  this  ratio.  No  cali¬ 
bration  is  performed  for  mean  spectral  period  or  any 
other  bulk  parameter  associated  with  spectral  distri¬ 
bution  of  energy. 


4.  Dissipation  model  on  parametric  spectra 

In  this  section,  to  further  examine  the  frequency- 
dependent  characteristics  of  the  calibrated  dissipation 
model,  we  analyze  calculations  with  the  new  observation- 
consistent  dissipation  term  on  parametric  spectra.  Whereas 
in  other  sections  the  time  evolution  of  spectra  is  modeled, 


computations  here  are  performed  outside  the  wave 
model,  and  interactions  with  other  source  terms  are  not 
considered. 


a.  Description  of  composite  parametric  spectra 

Because  the  directional  dependence  of  the  new  dis¬ 
sipation  function  (20)  only  exists  via  the  multiplication 
by  the  directional  spectrum,  that  is,  /Cds(/)  =  Sds(/,0)/ 
*he  computations  can  be  per¬ 
formed  on  omnidirectional  parametric  spectra.  There  is 
considerable  evidence  for  wave  frequency  spectra  ex¬ 
hibiting  two  subranges:  /“ 4  dependence  closer  to  spectra 
peak  and  f~s  dependence  at  higher  frequency  (Forristall 
1981;  Ewans  and  Kibblewhite  1990;  Hwang  and  Wang 
2001).  The  representative  parametric  spectra  is  modeled 
as  a  composite  function  of  both  slopes  separated  by  a 
transition  frequency,  which  is 


E(f)  = 


f^Don  </)• 

W/, r)  (j)  ^ 


f  —ftT 
f>ftx 


(25) 


Here,  our  composite  spectra  has  flT  =  3 fp  based  on  ref¬ 
erences  cited  above. 

The  £Don(/)  is  the  spectral  model  proposed  by 
Donelan  et  al.  (1985), 
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Fig.  5.  The  composite  parametric  spectra  described  in  the  text. 
All  models  have  the  same  wave  height  and  peak  period.  Vertical 
axis:  (f)f4.  Horizontal  axis:  f!fp. 


where  the  coefficients  a  Don*  YDon>  and  ^Don  are  the  equi¬ 
librium  range  (rear  face)  parameter,  the  peak  enhance¬ 
ment  factor,  and  the  peak  width  parameter,  respectively, 
and  can  be  determined  based  on  inverse  wave  age  f/10/Cp 
shown  in  Young  (1999).  This  is  very  similar  to  the  “Combi 
spectrum”  as  given  by  Tsagareli  (2009)  and  Babanin  et  al. 
(2010),  though  those  authors  use  /  =  2.5ghrUw ,  following 
Kahma  and  Calkoen  (1992).  The  composite  parametric 
spectra  are  illustrated  in  Fig.  5  for  wind  speed  Uiq  = 
12  m  s^1  and  f/10/Cp  of  0.9  and  3.5  corresponding  to, 
respectively,  fully  developed  and  young  seas. 

b.  Spectral  dissipation 

In  the  point  model  simulations,  the  spectral  distribu¬ 
tion  is  allowed  to  evolve  freely,  so  each  of  the  four 
models  produce  different  spectral  shapes:  this  is  evident 
from  the  different  peak  period,  mean  spectral  period, 
and  frequency  narrowness  in  Fig.  3.  However,  it  is  useful 
to  compare  computations  of  dissipation  term  on  iden¬ 
tical  spectra  to  contrast  the  unique  behavior  of  each 
model.  The  four  new  dissipation  functions,  applied  to 
the  composite  spectral  model  for  fully  developed  and 
young  seas,  are  shown  in  Fig.  6.  The  simulation  for  a 
fully  developed  sea  (Fig.  6a)  reveals  remarkable  differ¬ 
ences  between  the  dissipation  functions.  In  particular, 
the  £(/)  =  E(f)  dissipation  model  DL1M1  predicts  that 
most  dissipation  is  in  the  region  of  /  <  2/p,  whereas 


the  E(f)  =  ET(f)  model  UL4M4  predicts  that  most 
dissipation  will  occur  in  the  region  of  /  >  2/p.  We  further 
point  out  that  the  dominant  wave  breaking  is  zero  or  near 
zero  for  the  fully  developed  spectrum.  This  behavior  is 
consistent  with  what  occurs  in  the  real  ocean  (YB06). 

For  a  young  sea  (Fig.  6c),  all  four  functions  show  most 
dissipation  occurring  at  or  near  spectral  peak.  Of  the 
four  models,  the  UL4M4  model  demonstrates  a  signifi¬ 
cant  bimodal  distribution  in  frequency  space,  with  dis¬ 
sipation  peaking  near  the  spectral  peak  and  near  3 fp. 
Recall  that  all  models  were  calibrated  in  section  3  such 
that  the  frequency-integrated  7j  versus  T2  ratio  is  roughly 
consistent  between  them  at  the  quasi-equilibrium  state; 
the  ratios  of  T\  and  T2  seen  in  Figs.  6b, d  reflect  this.  This 
is  qualitatively  similar  to  the  maximum  breaking  rates 
reported  by  Babanin  et  al.  (2007b),  which  showed  an 
initial  strong  increase  with  frequency,  then  a  leveling 
off  or  slight  decrease.  However,  there  are  limits  to  such 
comparisons,  because  dissipation  rate  is  dependent  on 
not  only  breaking  rate  but  also  breaking  intensity. 

The  behavior  of  the  models  clearly  demonstrate  that 
it  would  be  a  mistake  to  conflate  T\  with  low-frequency 
dissipation  or  T2  with  high-frequency  dissipation:  the 
two  modes  coexist  for  all  frequencies  above  the  peak. 

5.  Model  behavior  in  free  simulations 

In  this  section,  we  utilize  the  same  duration-limited 
simulations  as  described  in  section  3.  Simulations  are  “free” 
insofar  as  spectra  are  time  evolving  and  not  prescribed  (cf. 
section  4  above).  The  behaviors  of  the  new  observation- 
consistent  input  and  dissipation  formulas  are  analyzed  here, 
first  by  evaluating  the  time  evolution  of  the  integrated 
source  terms,  and  then  comparing  duration-limited  growth 
curves  against  empirical  growth  curves. 

a.  Time  evolution  of  source  terms 

Frequency-integrated  source  terms  from  the  point 
model  simulations  are  shown  in  Fig.  7.  As  expected,  the 
DBYB  wind  term  is  stronger  for  younger  waves,  being 
calculated  from  ({VIC)  -  l)2  rather  than  ((UIC)  -  l)1, 
as  with  the  KHH  input  term.  After  the  period  of  initial 
growth,  all  of  the  new  models  show  a  slow  weakening 
of  the  wind-input  term  with  a  slow  strengthening  of 
the  dissipation  term.  The  former  is  due  to  the  effect  of 
the  downshifting  of  the  spectrum  on  wave  celerity  C. 
The  integral  of  the  nonlinear  terms  is  expected  near  zero, 
deviating  some  because  the  maximum  frequency  for  in¬ 
tegration  is  only  1.0  Hz.4  The  two  models  with  a  tendency 


4  Ideally,  the  computation  should  be  performed  using  an  accu¬ 
rate  nonlinear  solver,  for  example,  as  used  by  Tsagareli  et  al.  (2010), 
and  extended  into  high  frequencies  (e.g.,  10  Hz). 
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Fig.  6.  (a),(c)  Sds(/)  vs  flfp.  (b),(d)  Tt(f)  and  T2(f)  vs //f.  (a),(b)  Fully  developed  wave  conditions 
<yc,  =  0.9.  (c),(d)  Young  wave  conditions  U]0iCp  =  3.5.  The  models  shown  are  the  same  as  those 
in  Fig.  3.  Computations  are  on  the  parametric  spectrum  with  transition  at  ftT  ~3 fpy  and  are  for  the  case  of 
U\o  =  12  m  s~*. 


to  overpredict  high-frequency  energy,  DL1M1  and  KHH’s 
n  =  1,  both  show  especially  strong  negative  values  for 
!SJf)df;  apparently  the  nonlinear  term  is  partially 
compensating  for  weak  Sds  in  the  higher  frequencies  by 
transferring  energy  from  the  high-frequency  range  to 
the  diagnostic  tail.  All  of  the  new  models  show  similar 
values  for  D[o[  =  J  Sds(f)df  for  large  wave  age  values, 
but  during  the  initial  stages  of  strong  growth,  the  Dtol 
of  the  E(f)  -  ET(f)  models  are  larger  than  that  of  the 
KHH  n  =  2  model  (compensating  the  strong  input  from 
the  DBYB  5in),  whereas  that  of  the  £(/)  =  E(f)  model 
(green  dashed  line)  is  quite  weak  at  this  stage.  The  two 
Af  =  4  models  show  some  oscillation  of  at  these  early 
stages;  this  is  apparently  associated  with  high-frequency 
energy. 

The  lower  right  panel  of  Fig.  7  shows  the  ratio  be¬ 
tween  dissipation  and  input,  R  =  Dlol//lot  =  JSds(/)d// 
J  Sin(f)  df .  According  to  Donelan  (1998)  (see  also  Fig.  4 
of  Babanin  ct  al.  2010),  this  ratio  should  diminish  from 
0.85  at  Ul0fCp- 5.5  to  0.97  at  1.5  <  UJCp  <  4.5. 
Unlike  Babanin  ct  al.  (2010),  we  do  not  enforce  this 
ratio  as  a  physical  constraint  but  allow  it  to  develop 
freely.  (This  is  discussed  further  in  section  6.)  All  of  the 
new  models  predict  ratios  lower  than  that  suggested  by 


Donelan  (1998).  By  2  h  into  the  simulation,  all  models 
except  for  DL1M1  consistently  show  R  ~  0.7,  which  is 
still  lower  than  the  Donelan  (1998)  value  of  R  =  0.97. 
Increasing  the  frequency  range  of  computations  in¬ 
creases  the  ratio  only  modestly,  with  computations 
to  10.0  Hz  (not  shown),  R  =  0.92  after  12  h  with  the 
UL4M4  model,  compared  to  R  =  0.87  with  computa¬ 
tions  to  1.0  Hz  (shown). 

b.  Energy  growth  versus  empirical  expressions 


Traditionally,  empirical  growth  curves,  such  as  those 
of  Kahma  and  Calkoen  (1992)  nondimensionalized 
and  plotted  on  logarithmic  scaling,  are  used  to  calibrate 
wave  models  (e.g.,  Tolman  and  Chalikov  1996;  van  der 
Westhuyscn  et  al.  2007).  Such  growth  curves  were  not 
considered  in  our  calibration,  but  it  is  worthwhile  to 
compare  to  such  expressions  post  facto,  as  we  do  now. 
Duration-  rather  than  fetch-limited  comparisons  are  made, 
for  reasons  explained  in  section  3.  The  comparison  is 
shown  in  Fig.  8.  Here,  the  nondimensionalized  time 
is  (-gtlU.0  and  the  nondimensionalized  energy  is 
e  =  {Egt-fUi q),  where  E  is  the  total  variance  or  energy. 

As  in  Fig.  3,  we  include  the  KHH  n  =  1  and  KHH  n  = 
2  models  for  reference.  However,  rather  than  showing 
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FlG.  7.  Frequency-integrated  source  terms.  Models  shown  are  the  same  as  those  shown  in  Fig.  3. 
(bottom  right)  The  ratio  between  dissipation  and  input  vs  inverse  wave  age  is  shown,  with  circles  in¬ 
dicating  the  points  2  h  into  each  simulation.  The  frequency  range  for  the  integration  is  0.042-1.0  Hz. 


the  same  new  models  as  in  Fig.  3,  we  use  this  opportunity 
to  quantify  the  sensitivity  to  the  primary  coefficient  con¬ 
trolling  the  strength  of  nonbreaking  dissipation.  We 
select  a  single  model  (UL4M4)  and  vary  the  parameter 
Recall  that  fe  =  0  in  Fig.  3;  the  other  two  values  shown 
here  arc  fe  =  0.006  and  fr  =  0.011,  based  on  the  recom¬ 
mendation  by  Ardhuin  et  al.  (2009,  2010)  that  0.004  < 
fe  <  0.013  is  a  reasonable  range  of  values  in  the  large 
ocean  basins.  The  UL4M4  model  is  selected  for  this 
because,  based  on  comparison  in  Fig.  3  of  the  four  new 
models,  it  is  the  one  expected  to  have  the  fewest  problems 
with  overprediction  of  high-frequency  energy  in  subse¬ 
quent,  more  realistic  simulations.  All  numerical  model 
growth  curves  are  created  using  simulations  with  U\ o  = 
12  m  s~!. 

The  four  empirical  growth  curves  included  are  as 
follows: 

1)  A  linear  fit  to  the  tabulated  data  of  Stewart  (1961), 
e  =  8  X  10~9£12  within  the  range  10  X  103  <  £  < 
22.0  X  103. 

2)  The  Sanders  (1976)  growth  curve,  as  presented  by 
Young  (1999,  e.g.,  his  Fig.  5.14).  This  is  given  by  e  = 
3.22  X  1(T3  tanh2(1.26  X  10'3^0'75). 

3)  The  CERC  (1977)  growth  curve,  as  presented  by 
Young  (1999,  also  shown  in  his  Fig.  5.14).  Here,  the 


growth  curve  is  e  =  5.0  X  10“3  tanh2(0.0125*a42), 
and  Young  (1999)  converts  nondimensionalized  fetch 
X  to  nondimensionalized  time  £  according  to  an 
expression  from  CERC  (1977). 

4)  The  final  expression  given  by  Kahma  and  Calkoen 
(1992),  which  is  for  the  “composite  dataset’'  with  U\q 
normalization,  e  =  5.2  X  10"  y  9  We  convert  the 
expression  to  nondimensionalized  time  £  according 
to  the  same  expression  from  CERC  (1977),  as  men¬ 
tioned  above. 

Pierson  and  Moskowitz’s  (1964)  PM  limit  is  also  shown 
in  Fig.  8,  and  is  the  total  energy  calculated  from  in¬ 
tegration  of  the  well-known  PM  spectrum  given  in  that 
paper,  as  is  done  by  KHH,  which  gives  £pm  =  3.6  X  10“3. 
The  tabulated  values  of  Moskowitz  (1964)  are  also 
plotted,  though  it  should  be  noted  that,  according  to 
Moskowitz, 

the  durations  and  fetches  do  not  satisfy  the  theoretical 
requirement  that  they  describe  the  time  required  and  the 
distance  needed  to  generate  a  fully  developed  sea  start¬ 
ing  from  zero  wave  conditions.  For  the  higher  winds  in 
particular,  the  tabulated  values  often  represent  the  time 
and  the  distance  for  a  sea  raised  to  a  given  height  by  a 
wind  of  lesser  velocity  to  grow  to  full  development  at 
the  higher  velocity. 
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FlG.  8.  Verification  of  duration-limited  growth  rate  (energy  vs  time).  Both  quantities  are 
nondimensionalized  using  Ul0.  Two  models  based  on  KHH,  three  new  observation-consistent 
models  (UL4M4),  four  empirical  growth  curves  (described  in  the  text),  the  PM  limit  and  PM 
values  are  shown. 


- KHH  n=l 

- KHH  n=2 

- UL4M4,  fe=0.000 

- UL4M4,  fe=0.006 

-  -  -  UL4M4,  fe=0.011 
1 1  m  1 1 1  Stewart  (1961)  (fit) 

Mm  Sanders/Young 
nun  CERC/Young 

Kahma  and  Calkoen  (conv.) 
——PM  (1964) 

O  Moskowitz  (1964) 


In  other  words,  qualitatively,  the  plotted  points  should 
be  shifted  to  the  right  in  the  figure,  but  quantitative 
means  are  not  available.  A  similar  warning  could  be 
made  with  regard  to  the  Stewart  (1961)  linear  fit. 

An  interesting  outcome  of  the  comparison  is  that 
the  numerical  models,  including  the  KHH  models, 
match  each  other  much  more  closely  than  the  empirical 
curves  match  each  other.  In  some  sense,  this  vindi¬ 
cates  the  calibration  method,  because  had  the  Sanders 
(1976)/Young  (1999)  curve  or  Stewart  (1961)  linear  fit 
been  used,  the  calibration  may  have  been  unsuitable. 
On  the  other  hand,  the  numerical  models  match  the 
CERC  (1977)  and  Kahma  and  Calkoen  (1992)  curves 
rather  well. 

The  nonbreaking  dissipation  has  a  noticeable 
effect  on  energy  predictions  near  full  development. 
This  comparison  has  implications  for  the  treatment 
of  the  nonbreaking  dissipation  in  the  calibration: 
the  models  with  nonbreaking  dissipation  appear  to 
undershoot  the  Moskowitz  values.  This  suggests 
that  the  preliminary  calibration,  which  was  created 
using  fe  =  0,  may  lead  to  underprediction  of  energy 
near  full  development  in  subsequent,  realistic  simu¬ 
lation  with  nonzero  fe.  This  is  easily  addressed  by 
considering  alternate  calibrations  produced  using 
nonzero  fe. 


6.  Discussion 

In  section  5,  we  looked  at  the  ratio  R  =  Dlot  IIio{  = 
j Sin(f)df  as  predicted  by  the  single-point 
model,  interpreted  in  the  context  of  the  observations  as 
discussed  by  Donelan  (1998).  Given  Ntoi  =  J  Snl(/)  df  = 
0  for  any  model  of  the  entire  frequency  range,  the  wave 
growth  is  controlled  by  7(ot  -  Dtot,  which  permits  an  in¬ 
finite  number  of  possible  values  for  R  —  Dtot//lot.  How¬ 
ever,  our  implemented  constraint  on  J[(S)in(/)/C]  df  in 
Eq.  (12)  implies  a  soft  constraint  on  Itot,  with  some 
variation  allowed  because  of  the  distribution  of  energy 
between  slower  and  faster  waves.  (This  variation  re¬ 
sulting  from  C"1  can  only  be  minor  for  a  model  that 
accurately  predicts  Tm0 i,  Tp ,  and  frequency  width.)  The 
effective  constraint  on  /tot  implies  that  only  minor  var¬ 
iation  in  Dtot  is  allowed  if  the  model  is  to  predict  wave 
height  accurately  (e.g.,  upper  left  panel  in  Fig.  3).  There¬ 
fore,  /tot  and  Dlot  are  already  essentially  fixed  without 
applying  a  constraint  on  Ry  and  any  model  that  does 
have  a  constraint  on  R  would  be  overconstrained  and 
thus  would  almost  certainly  fail  to  find  a  self-consistent 
solution.  The  modest  discrepancy  between  the  model- 
computed  R  values  and  the  Donelan  (1998)  values  sug¬ 
gest  a  shortcoming  either  in  the  implemented  constraint 
(/tot)  or  on  the  verification  via  Donelan’s  Dtot//lot  or  in 
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the  model  calibration  (/tot  -  Dtot).  This  is  left  as  a  point 
for  future  inquiry. 

As  mentioned  in  section  2b,  directional  narrowness 
parameter  A  is  set  to  unity  in  dissipation  calculations 
(19),  (21),  and  (22).  Implementation  of  dependence  of 
Sds(f) on  directional  spreading  may  be  pursued  in  a  later 
study  by  using  A  ^  1.  However,  because  the  directional 
narrowness  A  will  tend  to  decrease  with  increasing  fre¬ 
quency  in  the  spectral  tail,  this  will  reduce  dissipation  in 
the  tail,  which  could  present  new  problems  in  the  con¬ 
text  of  the  DBYB  input  term  with  its  relatively  strong 
forcing  in  that  region.  Further,  with  regard  to  (19)  [the 
calculation  for  Er  (/)]  Babanin  (2009,  see  discussion  of 
Fig.  5.28)  has  shown  that  the  parameter  is  probably 
unnecessary  in  any  model  that  includes  cumulative  dis¬ 
sipation. 

The  expression  (6)  is  given  by  DBYB  with  U  =  f/10. 
However,  it  is  desired  that  the  model  should  scale  with 
friction  velocity  £/*  rather  than  £/10.  The  reader  is  re¬ 
ferred  to  the  concise  summary  of  sealing  arguments  by 
KHH  and  the  issue  is  discussed  in  greater  detail  by  Alves 
et  al.  (2003).  In  earlier  versions  of  our  implementation, 
we  used  U  =  28f/+  in  (6),  similar  to  the  approach  taken 
by  KHH  with  the  intention  of  forcing  the  model  to  scale 
with  f/+.  However,  it  was  determined  experimentally 
that  this  modification  has  insignificant  impact  on  model 
results.  The  insensitivity  is  presumably  due  to  the 
physical  constraint  on  total  stress,  which  is  in  terms  of 
f/+,  and  applied  after  (6).  Thus,  for  the  representative 
wind  speed,  the  present  implementation  as  described 
in  section  2  uses  U  =  f/10,  consistent  with  DBYB.  The 
model  nevertheless  scales  well  with  £/*,  as  mentioned  in 
section  3. 

As  noted  above,  the  DIA  method  is  used  to  estimate 
the  four-wave  nonlinear  interactions  Snl  in  simulations 
herein,  even  though  the  single-point  model  computa¬ 
tions  are  certainly  feasible  with  more  complete  approxi¬ 
mations,  such  as  those  used  by  Tsagareli  (2009),  BY05 
and  YB06.  This  was  done  to  provide  consistency  with  a 
forthcoming  manuscript,  which  will  involve  more  ex¬ 
pensive,  two-dimensional  applications.  Also,  we  draw  a 
distinction  between  foundation-building  emphasis  of 
those  earlier  works  versus  the  present  work  with  em¬ 
phasis  on  practical,  routine  application.  The  less  accu¬ 
rate  nonlinear  interaction  computations  in  the  present 
simulations,  combined  with  the  high-frequency  limit  of 
the  prognostic  spectrum  of  1.0  Hz,  implies  that  conclu¬ 
sions  about  the  detailed  behavior  of  the  model  physics 
should  be  regarded  with  suspicion  until  confirmed  to 
also  occur  with  more  exact  computations.  For  example, 
it  will  be  shown  in  a  subsequent  manuscript  that,  in  com¬ 
parison  to  observations,  bias  in  frequency  spreading 
may  be  a  common  result  from  using  the  DIA. 


The  nonbreaking  dissipation  mechanism  of  Ardhuin 
et  al.  (2009)  is  used  herein  to  represent  the  slow  dissi¬ 
pation  of  nonbreaking  waves,  a  situation  represented  in 
the  model  as  E(f)  <  Er  (/).  The  Ardhuin  et  al.  method 
is  based  on  momentum  losses  by  the  waves  to  the  at¬ 
mospheric  boundary  layer.  Another  theorized  mode  of 
energy  loss  is  from  turbulence  generated  by  the  non¬ 
breaking  orbital  motion  of  waves  (Babanin  2006).  If 
the  viscosity  of  water  is  not  zero,  then  this  turbulence 
must  exist  (Phillips  1961).  This  and  other  methods  for 
accounting  for  nonbreaking  dissipation  will  be  evalu¬ 
ated  and  contrasted  with  the  Ardhuin  et  al.  method  in 
a  separate  study. 

The  implications  of  the  qualitative  behavior  demon¬ 
strated  in  Fig.  2  have  been  investigated  thoroughly  in 
numerical  experiments  here  (e.g.,  Fig.  6),  but  this  be¬ 
havior  can  also  be  discussed  in  an  intuitive  context. 
When  the  ratio  E(f)IEr(f)  changes  from  5  to  10,  for 
example,  as  a  result  of  current  shear  or  by  massive  wind 
input,  how  should  the  model  respond?  The  concave 
down  model  suggests  that  at  E(f)IET(f)  -  5  the  dissi¬ 
pation  Sds(f)  is  already  near  an  asymptote  and  does  not 
respond  strongly  when  E(f)fEr(f)  increases.  The  con¬ 
cave  up  models,  on  the  other  hand,  suggest  an  explosive 
increase  in  Sds(/),  especially  for  L  -  4  (or  M  -  4).  This 
scenario  may  be  interpreted  as  being  more  physically 
plausible  in  situations  where  the  wave  energy  must  be 
destroyed  within  small  time/space  scales  by  some  means, 
for  example,  in  the  surf  zone,  or  in  wave  blocking  situa¬ 
tions.5 

The  cumulative  dissipation  term  implemented  in  (22) 
represents  only  one  of  two  apparent  mechanisms  for 
cumulative  dissipation.  Specifically,  it  represents  the 
breaking  of  relatively  short  waves  in  the  wake  of  or 
on  top  of  large-wave  breaking,  triggered  by  the  large 
breaker.  However,  it  is  recognized  that  there  also  exists 
a  straining  action,  modulation  of  shortwave  trains  by 
the  underlying  large  waves,  which  causes  the  shortwaves 
steepness  to  increase  at  the  forward  faces  of  longer  waves, 
resulting  in  their  frequent  breaking  (e.g.,  Unna  1941; 
Longuet-Higgins  and  Stewart  1960;  Phillips  and  Banner 
1974;  Donelan  et  al.  2010).  This  type  of  dissipation  is 
obviously  distinct  from  (22)  because  it  does  not  require 
that  the  longer  waves  are  breaking,  and  it  was  in  fact 
implemented  experimentally  in  a  previous  version  of 
SWAN,  the  so-called  Cumulative  Steepness  Method 
(CSM;  see  van  Vledder  and  Hurdle  2002;  Hurdle  and 


5  We  point  out,  however,  that  none  of  these  models  have  been 
extended  for  use  in  the  surf  zone,  where  it  may  be  necessary  to 
introduce  new  features  to  the  dissipation  formulation.  The  reader 
is  referred  to  recent  work  on  this  topic  by  Filipot  et  al.  (2010). 
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van  Vledder  2004),  where  the  dissipation  is  estimated 
from  an  integration  of  the  steepness  of  longer  waves 


W0)  =  4  fartm/v/'W.*).  (27) 

.Jo 


where  <p(/')  is  a  function  to  produce  a  dependency 
on  the  relative  angle  between  £(/,  6)  and  the  lower- 
frequency  energy  acting  on  it,  E(f\0').  (In  that  earlier 
study,  the  CSM  was  implemented  as  the  sole  deep-water 
dissipation  mechanism.)  Prior  observations  of  the  sup¬ 
pression  of  wind  sea  by  swell  has  primarily  been  for 
cases  of  steep  swells  in  the  laboratory  (e.g.,  Phillips  and 
Banner  1974),  while  in  contrast  it  is  difficult  to  detect  the 
suppression  in  the  open  ocean  (Violante-Carvalho  et  al. 
2004)  where  swells  are  typically  less  steep.  Further,  it  is 
possible  that  the  suppression,  when  it  does  occur,  may 
be  due  to  the  reduction  of  wind  stress  rather  than  en¬ 
hanced  breaking  (Chen  and  Belcher  2000),  in  which  case 
the  additional  dissipation  term  would  not  be  appropri¬ 
ate.  Therefore,  it  was  decided  that  this  third  breaking 
mechanism  should  not  be  included  in  the  present  nu¬ 
merical  model  until  observational  studies  provide  clearer 
guidance. 


7.  Summary  and  conclusions 

This  paper  introduces  wind-input  and  whitecapping 
source  terms  consistent  with  key  features  observed  ex¬ 
perimentally  by  Banner  et  al.  (2000),  DBYB,  BY05,  and 
YB06.  The  source  functions  were  initially  developed  by 
Tsagareli  (2009),  Tsagareli  et  al.  (2010),  and  Babanin 
et  al.  (2010)  and  are  adapted  and  improved  here  for 
practical  application,  and  implemented  in  an  experi¬ 
mental  version  of  the  SWAN  model  of  Booij  et  al.  (1999). 
The  model  development  strategy  is  to  rely  on  observa¬ 
tional  studies  as  much  as  is  practical,  as  opposed  to  using 
theoretical  works,  such  as  Miles  (1957)  theory  (section  1), 
as  a  starting  point.  Until  now,  these  new  observation- 
consistent  terms  with  a  number  of  significantly  new  phys¬ 
ical  features  had  not  been  tested  in  a  two-dimensional 
model,  and  thus  have  not  been  available  in  a  form  that 
is  usable  in  practical  applications. 

The  new  physics  model  is  calibrated  and  evaluated 
using  simple  methods.  The  wind-input  term  is  taken 
from  DBYB,  which  is  based  on  their  observational  work, 
and  modified  here  by  including  a  physical  constraint  on 
the  total  stress  comparable  to  that  of  Tsagareli  et  al. 
(2010).  Also,  a  drag  coefficient  formulation  based  on  re¬ 
cent  observational  work  is  implemented.  The  calibration 
of  the  dissipation  model  herein  is  very  limited,  primarily 
concerned  with  model-to-model  consistency  rather  than 
direct  calibration  against  any  particular  observational 
dataset.  However,  the  new  dissipation  function  is 


observation  consistent  insofar  as  it  conforms  to  two 
features  of  dissipation  in  the  real  ocean  reported  in  the 
literature  during  the  past  decade.  The  first  feature  is  that 
dissipation  is  two  phase,  with  waves  of  any  particular 
frequency  dissipating  due  to  either  1)  the  instability  (and 
breaking)  of  waves  of  that  frequency,  or  2)  the  destabi¬ 
lization  by  larger  breaking  waves  (e.g.,  through  tur¬ 
bulence).  The  second  feature  is  that  wave  breaking  is 
thresholded,  such  that  when  the  local  spectral  density 
falls  below  a  spectral  threshold  derived  from  measure¬ 
ments,  no  breaking  occurs  at  that  frequency.  Though  they 
have  been  implemented  separately  in  some  fashion  pre¬ 
viously  (section  1),  these  two  features  were  not  included 
together  in  any  numerical  model  until  very  recently  [e.g., 
Tsagareli  (2009)  in  an  academic  model  and  Ardhuin  et  al. 
(2010)  in  WAVEWATCH  III].  Both  features  contrast 
sharply  with  dissipation  terms  of  the  previous  genera¬ 
tion  (e.g.,  Hasselmann  1974;  KHH;  WAMDI  Group 
1988;  Booij  et  al.  1999)  for  which  all  waves  are  considered 
breaking  all  the  time,  and  with  all  wave  systems  affect¬ 
ing  the  strength  of  dissipation  of  all  other  systems  in  a 
physically  implausible  manner  (Rogers  et  al.  2003). 

During  the  selection  of  calibration  coefficients  (dj,  02) 
for  the  new  dissipation  function,  neither  the  frequency 
distribution  (e.g.,  peak  period,  mean  spectral  period, 
frequency  narrowness)  nor  the  directional  distribution 
(e.g.,  directional  width)  are  considered.  However,  the 
dissipation  formulation  of  Babanin  et  al.  (2010)  is 
generalized  here,  permitting  variation  in  the  degree  or 
manner  of  sensitivity  to  the  breaking  threshold  exceed¬ 
ence.  Four  variants  (DL1M1,  UL2M2,  UL1M4,  and 
UL4M4)  are  selected  for  evaluation  here. 

Though  the  primary  intent  of  this  paper  is  to  introduce 
the  new  source  functions  as  implemented  in  a  practical 
model,  the  following  conclusions  can  be  made  from  the 
simple  computations  herein: 

•  The  strong  input  to  high-frequency  waves  by  the 
new  wind-input  source  function  (DBYB)  necessi¬ 
tates  strong  dissipation  at  these  frequencies.  The 
£(/)  =  £(/)  dissipation  model,  denoted  DL1M1,  in¬ 
sufficiently  damps  these  waves. 

•  It  is  sufficient  to  treat  the  calibration  coefficients  of  the 
new  dissipation  terms  as  constant,  that  is,  independent 
of  wave  age.  None  necessitate  complex  calibrations,  such 
as  the  lookup  table  procedure  discussed  in  section  3. 

•  All  four  variants  of  the  new  model  exhibit  smaller 
ratios  of  integrated  dissipation  to  integrated  wind 
input  Dtot//tot  than  those  reported  by  Donelan 
(1998).  This  indicates  that  the  models  retain  a  relatively 
larger  portion  of  wind  energy  input  than  given  by 
Donelan  (1998).  This  discrepancy  is  being  investi¬ 
gated  in  detail  in  a  separate  study. 
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•  In  comparisons  that  allow  assessment  of  model  quality, 
the  three  E(f)  =  ET(f)  models — denoted  UL2M2, 
UL1M4,  and  UL4M4 — show  very  similar  behavior, 
implying  a  surprisingly  weak  sensitivity  to  L  and  M 
in  (21)  and  (22)  (Figs.  3  and  7). 

•  The  UL4M4  model  adequately  follows  empirical  tem¬ 
poral  growth  curves  (Fig.  8),  though  the  comparison 
also  suggests  that  calibration  may  need  to  be  updated 
to  weaken  breaking  dissipation  using  wave  observa¬ 
tions  and  realistic  simulations  of  large  ocean  basins  with 
the  nonbreaking  dissipation  included. 

•  The  four  models  exhibit  widely  different  distributions 
of  inherent  dissipation  T\  and  induced  dissipation  T2 
with  frequency,  even  when  computed  on  identical 
parametric  spectra.  In  such  a  comparison,  the  three 
E(f)  =  ET(f)  models  show  strong  sensitivity  to  L  and 
M  in  (21)  and  (22)  (Fig.  6).  However,  existing  observa¬ 
tions  do  not  provide  sufficient  means  to  assess  relative 
quality  of  these  predictions.  This  highlights  a  challenge 
to  address  in  further  research. 

In  a  forthcoming  manuscript,  wc  will  further  evaluate  the 
three  E(f)  =  ET(f)  models  using  a  variety  of  regional- 
scale  hindcasts:  cases  representing  pure  wind  sea,  mixed 
sea/swell,  hurricane,  and  slanting-fetch  conditions. 
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